0B 



MATHEMATICAL ASPECTS 
OF MODELS FOR 
HAMILTON HARBOUR, 1978 



HAM 
MAT 

ANM. 




The Honourable 
George R. McCague, 
MinistP/ Minister 

of the K.H.Sharpe, 

Environment Deputy Minister 

io 



Copyright Provisions and Restrictions on Copying: 

This Ontario Ministry of the Environment work is protected by Crown copyright 
(unless otherwise indicated), which is held by the Queen's Printer for Ontario. It 
may be reproduced for non-commercial purposes if credit is given and Crown 
copyright is acknowledged. 

It may not be reproduced, in all or in part, for any commercial purpose except 
under a licence from the Queen's Printer for Ontario. 

For information on reproducing Government of Ontario works, please contact 
ServiceOntario Publications at copyri ght@ontario .ca 



•i 



i^ 



A Developmental Study of Some Numerical and 

Mathematical Aspects of Models for 

the Hamilton Harbour 

by 



H. Rasmussen 

Department of Applied Mathematics 

University of Western Ontario 

London, Ontario 



April 1978 



This report describes the work carried out 
under Ontario Ministry of the Environment 
Purchase Order No. A50340. 



G<\f^V 



1. Introduction 

The effects of different technical solutions to water quality 
problems in lakes are often studied using a mathematical model of 
the lake. This model is then transformed into an equivalent 
numerical model which is usually solved on a large computer. 

There are two main problems associated with such numerical models. 
One is due to the fact that it Is not possible to describe exactly 
all the physical processes which take place in such complex 
systems. For example the turbulent processes are characterized by 
empirical parameters; the values of these must then be estimated 
from measured field data. In the first part of this report we 
discuss different methods for obtaining estimates of these para- 
meters; and present the conclusions that can be drawn from the 
application of them to the Hamilton Harbour Model, see the Hamilton 
Harbour Study (1974). More detailed treatments are given in 
Appendices A and B. 

The second problem is the amount of computing time and hence cost 
for solving these numerical models. In order to determine the 
feasibility of reducing cost we are examining some models for 
which at least part of the solution procedure can be done analyti- 
cally, so that the expensive numerical treatment would be reduced. 
In Appendix G such a model is derived. 

2. Model validation 

The output of a mathematical model is generally in the form of time 
series of water velocities and waste concentrations at some point 
of the lake. Such models contain unknown physical parameters, whose 
values must be adjusted so that the model output agrees reasonably 
well with the corresponding measured data. This process is called 
model validation or system parameter identification. 

It consists of two parts: a measure of how close the calculated 
data is to the measured data, and, secondly, a systematic procedure 



for adjusting the parameters so that the two time series become 
closer to each other. The reason that a mathematical definition 
of closeness of the two sets of data Is required, is that in 
practical situations they will never be identical. In this report 
we mainly use the square of the differences and maximum difference 
of the two sets of data as our measures of how close the calculated 
data is to the measured data. 

The main conclusion that can be drawn from the results presented 
In Appendix B is that with the given time step, 20 sec, and the 
given space mesh the major features of the flow in Hamilton Harbour 
are reasonably well modelled, while the more detailed features are 
not as well reproduced. However it seems likely that with a smaller 
time step and a finer mesh a better description of these detailed 
features would be produced at an increased computing cost. 

3. Semi-analytic methods for lake modelling 

One of the problems with the standard numerical models of lakes, 
such as those devised by Leenderste (1970), is the large amount of 
computing time required to solve them. It seems reasonable to 
expect that one should be able to reduce this cost by using models 
which can be partially solved analytically so that a computer would 
only be required for a smaller part of the solution procedure. 

The Ekman-type model which is based on the assiunption of small 
Rossby numbers and small horizontal turbulent mixing satisfies this 
criterion. It has been applied by for example Liggett and 
Hadjltheodorou (1969) and Gedney and Lick (1972). 

A careful study has been done of the derivation of this model, 
and in Appendix C it is shown that the standard formulation of the 
boundary conditions must be modified. This should lead to more 
accurate predictions by the model of the water velocity. We plan 
to apply a model of this type to the Hamilton Harbour. 
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Appendix A 

Review of some methods for 
validation of lake models. 

1. Introduction 

In the formulation of mathematical models of the circulation and 
transport of pollutants in lakes it is not possible to exactly 
describe all the physical processes which take place in such complex 
systems. For example at present the turbulent processes can only 
be described approximately by unknown parameters multiplying the 
gradients of quantities such as mean velocity components and mean 
pollutant concentrations. Thus one is faced with the problem of 
attaining values of these empirical parameters. Usually rough 
estimates can be obtained from laboratory experiments, but in order 
to get more exact values the model output must be adjusted until 
it compares reasonably well with extensive field data. A review of 
some of the procedure for solving this type of problem, which is 
often called system parameter identification, form the subject of 
this appendix. 

Only little work seems to have been reported in the literature on 
the systematic estimation of the empirical parameters in lake 
models. In a series of papers Simons (1974, 1975, 1976) considered 
the problem of verifying a numerical model of Lake Ontario. Since 
his model of the fluid flow was linear he could estimate an effective 
wind-stress coefficient by the comparison of observed and computed 
water levels and vertically mean currents. The bottom stress 
coefficient was taken to be equal to the wind-stress coefficient. 
The presented comparisons in Simons (1974) are for three different 
time scales: a) currents averaged over three days, b) current 
smoothed by a digital filter which removed all frequencies equal 
to or higher than the frequency of the inertial oscillation for Lake 
Ontario, and c) periods less than or equal to the inertial 
period. Only for the last time scale are the effects of varying 
some of the turbulent parameters presented. 



In Penumalli et al (1976) estimates of the turbulent dispersion 
coefficients were obtained for an estuary in Texas by adjusting 
these coefficients in a steady state model until the least square 
difference between the model output and the measured field data is 
sufficiently small. In this study the components of the water 
velocity vector were prescribed and the mathematical model consisted 
of the linear mass balance equation. 

In this report we describe some of the techniques developed in other 
fields of engineering for system parameters identification and 
unconstrained optimization and show how they can be applied to a 
typical model of a lake. These techniques have been developed and 
applied to chemical engineering, [Seinfeld (1969)], hydraulics, 
e.g. [Yeh (1975)] and other fields. A general review of the 
application of some of these techniques to control engineering is 
given in Astrora and Eykhoff (1971) . 

The importance of the verification of a given lake model lies in 
the fact that if the model is not correctly verified, at best only 
very qualitative results can be obtained. 

2. Formulation of a specific model 

In order to illustrate the discussion in this appendix we will consider 
the depth-integrated two-dimensional model for water movements used 
in the Hamilton Harbour Study (1974). 

Let (x,y) be a cartesian coordinate system with t being the time 
variable. Then the model consists of the equations of motion 

2 2 k X 

U, + UU + VU - fV + g c + g "^^ t ^ ^ -TT = (1) 

t X y C H ^ 

2 2 ^ V 
Vj. + UV^ + W^ + fU + g Cy + g ^^" ^ ^ ^' - ^ = ^ ^2> 

C H 



and the continuity equation 

HTT\ + 

9y 



^t "^ ^ (»u> -^ -k ^"^^ = ° <3> 



where U = 3U/Sx, etc. 

X 



U = velocity in x direction 

V = velocity in y direction 

f = Coriolis parameter 

g = acceleration due to gravity 

C = Chezy coefficient representing bottom friction 

H = depth = C + reference depth 

c, = elevation of water surface above reference 

T = component of the wind-stress in the x direction 

y 

T = component of the wind-stress in the y direction 
p = density of water. 

X y 

The relationships between the wind-stress terms t and t and the 

wind are 



T = e w sin ^^^ 

y 2 
T = w cos ij; 



where 9 = wind-stress coefficient 
w =» wind speed 
i|j = direction of the wind. 

This model contains two unknown parameters: the Chezy coefficient 
C representing bottom friction and the wind-stress coefficient 8. 
Some systematic method for choosing these parameters must be used 
if the output of the model is to be trusted. In this paper we 
consider minimization of a criterion function. In the remaining 
part of the paper it will be assumed that a numerical procedure for 
solving equations (1), (2) and (3) together with appropriate boundary 
and initial conditions has been developed, 

so that accurate numerical approximations to U, V and C can be 
calculated. 

3. Minimization of a criterion function 

In this method the unknown parameters C and 6 are chosen such that 
a measured variable, for example the velocity at a point, agrees 



well with the calculated values of the same quantity. This raises 
three questions: 

1) What quantities should be used for this comparison? 

2) What criterion should be used to measure how well 
observed and calculated quantities agree? 

3) How can a systematic procedure be devised for adjusting 
C and 9 so that the criterion in 2) is satisfied? 

We shall now consider these three points in turn. 

1) Choice of comparison quantities. 

There are basically four different types of quantities which can be 
used for comparison purposes (see Seinfeld (1969)). 

a) Time and space dependent quantities, e.g. the elevation c; at 
several points and at different times. 

b) Time dependent quantities, e.g. the volume by which the lake 
differs from the initial volume or the velocity at a given point. 

c) Space dependent quantities, e.g. the maximum elevation for the 
time interval of interest at each point of the lake. 

d) Time and space independent quantities, e.g. the maximum eleva- 
tion for the time interval of interest for all points of the lake. 

If the measured data is already available, the choice of comparison 
quantities is made. However if the data has not been collected at 
the start of the program, careful thought should be given to the 
choice of comparison quantities within the technical limitations of 
the equipment available for measurement. In general one is restricted 
to velocity measurements at different points at regular time 
intervals. It is important to obtain measurement for reasonably 
long time period. 

2) Criterion functions 

The two most commonly used criterion functions are 

B 



a) the least square criterion, and 

b) the minimax criterion. 

a) Least square criterion 

Suppose the measured data consists of the velocities U and V at 
a fixed point (x*, y*) at regular time intervals t., j = 0, 1, . 
N; write these measurements as 



Uj^(tj) = UCx*, y*, tj) 



(5) 



Vj^(tj) = V(x*, y*, tj) 



The corresponding calculated quantities we denote by 

U^(t , C, e) and V^(t , C, 9); 

these, of course, depend on the particular values of C and 9 used 
in the calculations. We now define 



D(c,e) = I {[u^(tj) - u^(tj,c,e)]' 



+ [V^(t.) - V„(t,,C,9)]^} (6) 

The problem is now to find those values of C and 9 for which D is 
a minimum, i.e. 



min D(C,9) (7) 

c,e 

This is the classical least squares criterion with uniform weights, 
b) Minimax criterion 
If we define 



G, = lu (t ) - U^(t ,c,e) 



"^^ = V^^^i^ - V_(t ,C,9) 
J m j C j 



(8) 



the minimax criterion function can be written in the form 



min max (c + n,); j = 1, 2, ..., N 

c.e j ^ J 

here we have assumed uniform weighting functions. This can also 
be written in the form 



min D(C,e) (9) 

c,e 

where D(C,6) = max (e. + n.) 

j 3 J 

J = 1, 2. ..., N 

3) Numerical methods for calculating C and 9. 

The problem of calculating C and 6 is a typical unconstrained opti- 
mization problem. Many different procedures have been devised and 
presented in the literature. The following discussion is mainly 
based on Kowalik and Osborne (1968) . 

We shall consider three different methods and try to point out 
advantages and disadvantages of each method. One point which must 
be kept in mind is that usually it is very time-consuming to calcu- 
late D(C,6), so it is necessary to restrict our attention to methods 
which converge reasonably quickly. 

(i) The method of steepest descent. 

(k) (k) 
Suppose we have estimates for C and 9, say C and 9 , and we 

wish to improve on these estimates. We can do this by first 

(k) (k) 
computing the gradient of D at C and 6 . Numerically this 

gradient can be approximated by 

,(k) ,(k), „,^(k-l) Jk) 



C(k) ,(k) _ ^.(k-1) 

e 

(k) -(k), ^,^(k) ^(k-i) 



(k) ^ 

(k) 

: D(C^"\e^"0 - D(C^"\9^"""^) 
where i and j are mutually perpendicular unit vectors. 
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Thus we can write the new values of C and 9 as 



c(k+i) . c^k) ^ ,, . ,0 



e^'^-^^^ =e^^^ +aj . VD 



where a is a positive parameter. The question is now how to 
choose a. One way is to pick a so that 



,(,(k+l)^3(k+l)j 



is a minimum. It is also possible to give a a constant value at the 
beginning of the procedure and then use this value for each 
iteration. 

The advantage of the method of steepest descent is its simplicity; 
it is easy to understand and to apply. However, Its performance is 
many times disappointing; after a few satisfactory iterations the 
convergence rate become very small. It is not much used anymore. 

(ii) The method of Davidon, Fletcher and Powell. 

This method, which is also called a variable metric method, seems 
to be the best general-purpose optimization procedure making use 
of derivatives that is currently available. It is a modified 
gradient method; it uses certain information which is generated at 
each iteration to construct the Hessian matrix of the function. 
This indicates that the method attempts to use secondnarder infor- 
mation while the method of steepest descent uses only first-order 
information. 

In terms of a general function F(x) where x is the vector 

(x, , x», ..., X ) the i-th stage of the algorithm can be expressed 
12 n 

as follows. 

Compute _d = -H ^ 

compute ^. to minimze F(x + ^d.) 



1 



-i' 



(i+1) (1) ^ . . 
set X = X + ^j^ Ij^ 
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compute ^ = VF(x^ ) 



set 



T 



d d 

compute H, = H, , + ;^ ^ ~^ 



i-1 i T „ 

% h-l % 

T 
Here d Is the transpose of d^. The first stage of the algorithm 

(0) 
consists of choosing a starting point k and an initial approxi- 
mation H for the Hessian matrix; it is customary to set 



«0 = ^ 
where I is the unit matrix. It is also necessary to calculate 

For the particular problem we are concerned with n - 2 and 
x^ = C, x. = 6 and 

F(x) = D(C,e). 

The advantage of this method is its fast convergence rate, consi- 
derably faster than the method of steepest descent. However, it 
does require the calculation of the gradient of D or an approxi- 
mation to it which can be time-consuming. 

(ili) Powell's method without derivatives. 

In contrast to the two methods discussed above this method does 
not require the calculation of the gradient of D. It is based on 
the idea of conjugate directions. Again we will present the method 
for a general function F(x) where x is the vector (x^ , x^ , ..., 
X ). We assume that we initially have a starting point x and n 
independent vectors d^^, d2, .... d . Then the ith stage is 

Compute ^r> to minimze F(x + Ad ) 
U — — n 
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— — -T 



For j = 1, 2, ..., n compute A, to minimize 

F(x^^^ + Xd.) 

and set x^^"^""-^ = x^^^ + X. d. 

- - J -J 



set d = d.^j^. j = 1. 2, .... n-1 

(n+1) (n) 
set d = X - X 



(0) (n+1) 
set X = X 



Go back to the beginning. 

This method usually requires more iterative steps before convergence 
Is attained than the method of Davldon, Fletcher and Powell. However 
since it is not necessary to evaluate the gradient of D at each 
step, the total computing time required is often similar. 

For all three methods it is necessary to have a criterion for 
terminating the calculations. Usually it is that 

j,^c(n+l)^3(n+l)j -D(C^">.e(">) <. 

where £ is some preasslgned tolerance. A reasonable value is 
often 

-3 

G = 10 . 

4. Conclusions 

The desired sequence of steps in producing a reliable numerical 
model can now be stated as follows: 

a) sufficient field data is collected so that the main physical 
processes in the lake can be stated. 
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b) an adequate mathematical model containing various empirical 
parameters is formulated. 

c) a numerical procedure for solving this model is devised. 

d) field data is collected so that adequate boundary and initial 
conditions for the model are obtained. 

e) numerical experiments are carried out to ensure that the time 
and space step sizes are sufficiently small. 

f) adequate field data which is compatible with the model, is 
collected , 

g) a procedure, for example one of the procedures described in 
this report, is used to calculate reliable values of the different 
empirical parameters in the model. 

It is now possible to use the model with a fair amount of confi- 
dence to predict quantitatively the effects of several different 
technical solutions to water quality management problems. However, 
without the validation of the model being carried out usually only 
qualitative results can be expected. 
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Appendix B 

Application to the Hamilton Harbour model of 
some methods for the validation of lake models. 

1. Introduction 

Some of the procedures for model validation outlined in Appendix A 
were applied to the results obtained using the numerical model of 
Hamilton Harbour developed by the Water Resources Branch, Ontario 
Ministry of the Environment; see the Hamilton Harbour Study (1974). 
The results and the conclusions that can be drawn from them are 
presented in this appendix. The measured data used in these calcu- 
lations were supplied by the Water Resources Branch. The measured 
and calculated data were obtained for grid point (15,6). 

2. Formulation 

The model being considered in this appendix Is essentially the 
same as the model described In Appendix A except that the Chezy 
coefficient C which models bottom friction is replaced by 

C=i:^H^. 
n 

The parameter n, called Manning's n, can be assigned different 
values in different parts of the lake; in this study four different 
values of n, denoted by n^, n„, n_ and n. , were used. The numerical 
solution of the model was obtained using a program developed by the 
Water Resources Branch, Ontario Ministry of the Environment. 

3. Results 

Several solutions were obtained corresponding to different values 
of the wind stress parameter 9 and different sets of values for n, , 
n^, n~ and n, , Manning's n. Since the n's were changed by the same 
percentage only the value of n, is given for each solution; the 
base values are n, = .016, n^ = .023, n^ = .030 and n, = .040. 
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In Table 1 we present some results for 12 hours of real time using 

a time step of 20 sec. It is seen that the lowest value of the 

least square criterion is obtained at Runs 5 and 8 and the highest value 

at Run 15. Some of these results are also presented in Figure 1. 

This seems to indicate that the minimum of the least square criterion 

is in the neighbourhood of 6 = .0032 and n^^ = .0126 and that this 

minimum is not much less than 0.2374. Without extensive and very 

expensive computation it is not possible to state conclusively 

that this minimum is the smallest that can be obtained; i.e. it 

would be very expensive to prove that this minimum is global and 

not local. 

Some results for the max-min criterion are also given in Table 1 
and we see that the minimum is in the neighbourhood of 6 = .0035 
and n = .0126 with a minimum less than .0842. Again it would 
be very expensive to prove that this minimum is global and not 
local. 

Finally In Table 2 we give the mean values of the two velocity 
components. We see that only in Run 1 Is one of the velocity 
components in the wrong direction. 

From Table 1 we see that the optimum values for 9 and n^ are the 
ones used in Runs 5, 8 and 11. In order to analyse these runs in 
more detail we calculated the percentage differences between the 
mean and the root mean square of the measured velocities and the 
model output; these are presented in Table 3. 

4. Conclusions 

Several conclusions can be drawn from the results presented in the 
previous section. 

a) Major features seem to be modelled reasonably well. This is 
indicated by the fact that 

(i) For all runs except the first one the mean velocities are 
in the right directions. 
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(ii) The differences in the rms values of the measured and 
calculated data are not too large. 

b) With the present time step size and space mesh, the detailed 
features are not modelled as accurately as the major features; 
this conclusion is based on 

(1) The least square and min-max criteria are only satisfied 
very approximately. 
(ii) The difference in the rms values of the measured and 
calculated V velocities is around 50%. 

c) The best values for 9 and n. seem to lie in the neighbourhood 
of .0035 and .0126, respectively. However, Table 1 indicates 
that the two criteria give slightly different optimum values. 

It has been supposed the windstress parameter 6 does not vary from 
grid point to grid point but is constant. It is known that this 
is usually not the case, but since very little is known about this 
spatial variation of 9, it is difficult to incorporate this into 
the mathematical model. However the fact that it is not done may 
be one of the reasons why the detailed features of the flow in the 
harbour are not as well modelled as the major features. 
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Table 1: 


Least square 


and max crit 


Run if 


Wlndstress 


Manning's 




coefficient 


^1 


1 


.0020 


.012544 


2 


.0025 


.012567 


3 


.0030 


.016000 


4 


.0030 


.012609 


5 


.0032 


.012600 


1 


.0032 


.011970 


7 


.0032 


.012240 


8 


.0032 


.013000 


9 


.0032 


.013500 


10 


.0034 


.012591 


11 


.0035 


.012600 


il 


.0040 


.013400 


m 


.0040 


.011970 


m 


.0040 


.012559 


15 


.0050 


.012508 



I[(u -uj' 



m 



+ (v„ - V^) ] 
m c 



.3614 

.2951 

.2575 

.2447 

.2374 

.2410 

.2388 

.2374 

.2386 

.2493 

.2638 

.3387 

.4217 

.3844 

.5360 



Max[|u^ - U I 
m c 

'in r* ' 



.1034 
.0991 
.0952 
.0932 
.0898 
.0892 
.0894 
.0901 
.0906 
. 0861 
.0842 
.1043 
.1282 
.1205 
.1349 
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Table 2: Mean values of u and v for a 12 hour period. 



Run # 



Mean value 

of U 



Mean value 
of V 



I 

2 

3 

4 

5 

6 

7 

S 

9 
10 
11 
12 
13 
14 
15 
Measured 
data 



.01150 
.01601 
.02191 
.02193 
.02329 
.02277 
.02302 
.02353 
.02378 
.02380 
.02376 
.02421 
.01855 
.02133 
.01620 

.00712 



.00636 
-.00495 
-.01147 
-.01430 
-.01574 
-.01463 
-.01517 
-.01617 
-.01643 
-.01554 
-.01499 
-.01331 
-.00864 
-.01043 
-.00366 

-.02199 
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Table 3: Comparison of Runs 5, 8 and 11. 

Run 5 8 11 

B .0032 .0032 .0035 

nj^ .0126 .0130 .0126 

Least square 

criterion .2374 .2374 .2638 

Mln-max 

criterion .0898 .0901 .0842 

Difference in 

mean for U 227% 230% 234% 

Difference in 

rms for U 54% 54% 35% 

Difference in 

mean for V 28% 26% 32% 

Difference in 

rms for V 9% 10% 22% 
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Appendix C 

On the boundary conditions for 
wind-driven lake circulation models.* 

1. Introduction 

Approximate solutions for the steady wind-driven circulation in 
shallow lakes have been calculated by Hadjitheodorou (1) and 
Liggett and Hadjitheodorou (2) using an Ekman-type model based on 
the assumption of small Rossby number and small horizontal turbulent 
mixing. Later this analysis was extended by Gedney and Lick (3). 
A review of this type of model is given in Cheng, Powell and 
Dillon (4), 

The model can be expressed as a second order linear elliptic partial 
differential equation for a volumetric-transport stream-function ^. 
In (2), (3) and (4) the boundary condition imposed on t^r is that ti» 
is constant on all boundaries. In this paper we show that in 
addition to this condition, the boundary conditions on the velocity 
components imply that the normal derivative of ij^ must also vanish 
on the boundaries. However, for a second order elliptic equation 
both of these conditions cannot be satisfied simultaneously. It 
is shown in the last section of the paper how the model can be 
modified so as to remove this difficulty. From this we can draw 
the conclusion that this type of model is not adequate near the 
shore of a lake. In this region a more complete model containing 
all the viscous terras must be used. 

2. Formulation 

Let (x,y,z) and (u,v,w) be the nondimensional position and velocity 
vectors, respectively. In a rectangular coordinate system with x 
positive eastward, y positive northward, and z positive upward and 
zero at the surface. If p is the nondimensional pressure, the 



* This appendix has been written together with Mr. P. Forsyth, Jr., 
Department of Applied Mathematics, University of Western Ontario. 
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linearized equations for the wind-driven lake circulation are 

1 

-V = -p + TT U (1) 

X 2 zz 

2m 



u = -p + — =■ V (2) 

y 2m^ " 

P^ = -1 (3) 

u + V + w = (4) 

X y z 

2 
2 fD 
where m = -r — 

2n 

f = Coriollls parameter (assumed to be constant) 

D = typical vertical dimension 

n = eddy viscosity 

p = 9p/9x, etc. 

These equations are derived in Liggett and Hadjitheodorou (2). 



The boundary conditions are 

fLT fLT 



u = 



— = A, V = ^ = r, w = at z = (5) 



z ng z ng 

u = v»=w=0 at z = -h(x,y) (6) 



where 

z = -h(x,y) is the bottom of the lake 
L = typical horizontal dimension 
T = wind stress in the x direction 

X 

T = wind stress in the y direction 

y 

g = gravitational attraction 
3. Analysis 

Equation (3) implies that p and p are Independent of z, so 

X y 

equations (1) and (2) can be integrated to give 

13 I 



, ._ mz -mz. 
u = -p + cos niz(C_e - C.e ) 
y 2 4 

. ,„ mz „ — niz. ,_. 

- sin mz(C e - C-e ) (7) 

, ,_, mz , „ -mz, 

V = p + cos mz(C^e + C„e ) 

+ sin mz(C„e + C.e ) (8) 

where C, , C_, C- and C, are functions of x and y which are evaluated 

1 z J A 

by imposing the boundary conditions (5) and (6) ; they are given in 
the appendix. The pressure p is still unknown. 



An additional equation is obtained by the integration of the 
continuity equation (4) with respect to z from z = -h to z = 0; 
hence 

fS' 

u dz + V dz = (9) 



-h ^ J-h y 



which can also be written in the form 



^ udz + Y- vdz = (10) 

3^ J,h ^^ J-h 



since u=v=w=Oatz=-h and w = at z - 0. This equation 
is satisfied by the function ij;(x,y) defined by 



,0 

^ = - vdz and if/ 



X 



-r 



y j-h 



udz. (11) 



If we substitute u and v from equations (7) and (8) into the above 
expressions and integrate, we get 



ip = hh.p - hh,p - hh.,r + hh-A (12) 

X 4 X 1 y 3 2 

If; = hh,p + hh, p + hh„r + hh_A (13) 

y 1 X 4 y 2 3 

where h, , h„, h. and h, are functions of x and y given in the 

1 2 J 4 

appendix. These expressions can be solved for p and p , and since 

P = P (14) 

•^xy yx 
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the following equation for ^ is obtained 

y^^ = T7-(hJ + hf)[(r + s )4- + (r + a )i(; 
h^ 1 4 X y ^x y x y 

+ ^(qr + tA) + -^(tr - qA)] (15) 

where q, r, s and t are functions of h, h^^, h„, h- and h- given in 
the appendix. 



The analysis so far is similar to the one presented in Liggett and 
Hadjltheodorou (2). However the treatment of the boundary conditions 
on 'p is different. In order to simplify the subsequent analysis 
we shall consider a rectangular lake with boundaries parallel to 
either the x or y axes. Consider now one of the boundaries which 
we can assume to be along y = 0. As y approaches zero, we see 
from equation (11) that ip and ij) must also approach zero, since u 
and V must be finite, i.e. 

\b =1/) =Oaty = (16) 

X y 

These conditions, which also hold at the other boundaries, are 
equivalent to 

4j = C, lii =Oaty = (17) 

n 

where ]h denotes the normal derivative of i() and C is a constant. 

n 

If there is no inflow into the lake, we can take C to be zero. 
Liggett and Hadjltheodorou (2) and Gedney and Lick (3) only used 
the first of these conditions. 

Since equation (15) for ij; is a second order elliptic partial 
differential equation it is clear that we can only impose one of 
the two boundary conditions in (17) . It is easy to see how this 
problem of too many boundary conditions arose. In the derivation 
of the linearized equations (1) and (2) it was assumed that the 
horizontal viscous terms were much smaller than the vertical viscous 
terms and hence could be ignored. This effectively lowered the 
order of the equation; however, the complete set of boundary 
conditions (5) and (6) were retained. 



4. Treatment of the boundary conditions. 

In Liggett and Hadjitheodorou (2) and Gedney and Lick (3) it was 

implicitly assumed that the condition it = on the boundaries 

n 

could be ignored. However, as we can see from the definition of 
^, equations (11), this may lead to u and v being infinite at the 
boundaries. A more suitable treatment of the boundary conditions 
consists of supposing that the model defined by equations (15) and 
(17) is only valid in the interior of the lake. Along the shore 
there is a region in which a more complicated model is required. 

Let us now illustrate how this idea can be used in a numerical 
procedure. Suppose we have a rectangular lake with sides parallel 
to the X and y axes; say x=0, x=X, y=0 and y = Y where X and 
Y are given constants. We cover the lake by a mesh defined by 



x^ = i Ax, Yj = J Ay (18) 

where i = 0, 1, 2, .... X/Ax, and j = 0, 1, 2, .... Y/Ay. Consider 
now the region near y = 0, see Figure 1. At y = i(j = C and 
i<( = 0. The last condition can be approximated by 

to order Ay. Hence 

'^i.l = h,0 = ^' ^20) 

similar conditions can be obtained near the other boundaries. We 
have now a properly posed problem for the region 

1< i<X,i,i <j <JL_i (21) 

The difference between this problem and the problem solved by 
Liggett and Hadjitheodorou (2) and Gedney and Lick (3) is that 
here the depth function -h is not zero at the point where i> is 
equal to a constant. Since the grid spacings used in most finite 
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difference treatments of lake models are fairly large, e.g. Gedney 
and Lick (3) used x = y = 0.5 miles in some parts and x = y = 2 
miles in other parts, this difference can lead to large differences 
in the calculated current patterns, especially near the shore. 
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Appendix 

It is convenient to define several intermediate expressions due 
to the length of the equations. First, the following are defined: 

~ _ 2 , , -mh , mh. 2 , .2 , . -mh mh, 2 ,^„. 

a - cos mh(e + e ) + sin mh(e - e ) (22) 

$ = — ^ (sif^ "'h ~ cos mh) (23) 

J i_ / ~n>h mh. , „ , , 

Y = sin mh(e - e ) (24) 

1 — mH 

6 = -r- e (sin rah + cos mh) (25) 

zm ' 

, , -rah , rahv , „ , , 

e = cos mh(e + e ) (26) 

K = — e (sin mh + cos mh) (27) 

X = — e (sin mh - cos mh) (28) 

Then the expressions for C., C,, C and C can be written in the 
form 



C = c +-^ - A 
1 3 2m 2m 

^2 " " ^4 "^ 2^ "^ 2^ 

C3 = ^[-e|j+Y|^+ r(Sc - &y) + A(By + 5e) ] (29) 

C^-^[-Y|^-e|^+ r(BY + (5e) - A(Be + &y)] (30) 

The coefficients in equations (12) and (13) are 

h] = ;^ [y(6 + k) + e(6 + A)] (31) 



1 ah 



a(i + e - 6) 
1 r ni 



^2 = -^ [ 2^^ (By + 6e)(e + k) - (6e - 6y)(S + \)] (32) 

^3 " ^ [ "^^2m ^^ ^ ^^"^ ~ ^^^^^ + <) - (By + 6c) (6 + X)] (33) 
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h, = -r [- ah + e(B + <) - y(6 + \)] 
4 an 



(34) 



The coefficients In equation (15) are 



r = 



h(hj + hj) 



8 = 



h(h^ + h^) 



q = 



2 2 

h, + h, 

1 4 



t = 



^3\ ~ \^2 

2 2 



(35) 



Lake 



y = 



1-1 



i + 1 



J = 2 



j = 1 



j - 



shore 



Figure 1. Finite difference mesh near the shore. 
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